Black hole puncture initial data with realistic gravitational wave content 
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We present improved post-Newtonian-inspired initial data for non-spinning black-hole binaries, 
suitable for numerical evolution with punctures. We revisit the work of Tichy et al. [W. Tichy, B. 
Briigmann, M. Campanelli, and P. Diener, Phys. Rev. D 67, 064008 (2003)], explicitly calculating 
the remaining integral terms. These terms improve accuracy in the far zone and, for the first 
time, include realistic gravitational waves in the initial data. We investigate the behavior of these 
data both at the center of mass and in the far zone, demonstrating agreement of the transverse- 
traceless parts of the new metric with quadrupole-approximation waveforms. These data can be 
used for numerical evolutions, enabling a direct connection between the merger waveforms and the 
post-Newtonian inspiral waveforms. 
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I. INTRODUCTION 

Post-Newtonian (FN) methods have played a funda- 
mental role in our understanding of the astrophysical im- 
plications of Einstein's theory of general relativity. Most 
importantly, they have been used to confirm that the ra- 
diation of gravitational waves accounts for energy loss in 
known binary pulsar configurations. They have also been 
used to create templates for the gravitational waves emit- 
ted from compact binaries which might be detected by 
ground-based gravitational wave observatories, such as 
LIGO [U, and the NASA/ESA planned space-based 
mission, LISA I, g. However, FN methods have not 
been extensively used to provide initial data for binary 
evolution in numerical relativity, nor, until recently (see 
0, 0]), have they been extensively studied so that their 
limitations could be well identified and the results of nu- 
merical relativity independently confirmed. 

Until the end of 2004, the field of numerical relativ- 
ity had been struggling to compute even a single or- 
bit for a black- hole binary (BHB). Although debate oc- 
curred on the advantages of one type of initial data over 
another, the primary focus within the numerical rela- 
tivity community was on code refinement which would 
lead to more stable evolution. Astrophysical realism was 
very much a secondary issue. However, this situation 
has radically changed in the last few years with the in- 
troduction of two essentially independent, but equally 
successful techniques: the generalized harmonic gauge 
(GHG) method developed by Fretorius and the "mov- 
ing puncture" approach, independently developed by the 
UTB and NASA Goddard groups ^M. Originally in- 
troduced by Brandt & Briigmann [lO| in the context of 
initial data, the puncture method explicitly factored out 



the singular part of the metric. When used in numerical 
evolution in which the punctures remained fixed on the 
numerical grid, it resulted in distortions of the coordi- 
nate system and instabilities in the Baumgarte-Shapiro- 
Shibata-Nakamura (BSSN) [ll|, [13 evolution scheme. 
The revolutionary idea behind the moving puncture ap- 
proach was precisely, not to factor out the singular part 
of the metric, but rather evolve it together with the reg- 
ular part, allowing the punctures to move freely across 
the grid with a suitable choice of the gauge. 

A golden age for numerical relativity is now emerging, 
in which multiple groups are using different computer 
codes to evolve BHBs for several orbits before plunge and 
merger [H Q [H [H, [13, H H [13, IH ■ Comparison of 
the numerical results obtained from these various codes 
has taken place [13, [H, [13| , and comparison with FN 
inspiral waveforms has also been carried out with encour- 
aging success [I, [^ [1^ [1^ . The application of successful 
numerical relativity tools to study some important as- 
trophysical properties (e.g. precession, recoil, spin-orbit 
coupling, elliptical orbits, etc) of spinning and/or un- 
equal mass-black hole systems is currently producing ex- 
tremely interestingnew results [27 . 28, 29, 30,, S^, 
[H, [H, [H, [li, [13, M [H, [3, [U, [3 . It now seems that 
the primary obstacle to further progress is simply one of 
computing power. In this new situation, it is perhaps 
time to return to the question of what initial data will 
best describe an astrophysical BHB. 

To date, the best-motivated description of pre-merger 
BHBs has been supplied by FN methods. We might ex- 
pect, then, that a FN-based approach would give us the 
most astrophysically correct initial data from which to 
run full numerical simulations. In practice, FN results 
are frequently obtained in a form ill-adapted to numcri- 
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cal evolution. PN analysis often deals with the full four- 
metric, in harmonic coordinates; numerical evolutions 
frequently use ADM- type coordinates, with a canonical 
decomposition of the four-metric into a spatial metric 
and extrinsic curvature. 

Fortunately, many PN results have been translated 
into the language of ADM by Ohta, Damour, Schafer and 
collaborators. Explicit results for 2.5PN BHB data in the 
near zone were given by Schafer [4^ and Jaranowski & 
Schafer (JS) |44l |. and these were implemented numer- 
ically by Tichy et al. [i^. Their insight was that the 
ADM-transverse-traceless (TT) gauge used by Schafer 
was well-adapted to a puncture approach. To facilitate 
comparison with this earlier work we continue to 
use the results of Schafer and co-workers, anticipating 
that higher-order PN results should eventually become 
available in a useful form. 

The initial data provided previously by Tichy et 
al. already include PN information. They are accurate 
up to order {v/cf' in the near zone (r ^ A), but the 
accuracy drops to order (v/c)^ in the far zone (r 3> A) 
[here A ^ ny^ rf2/G{mi + TO2) is the gravitational wave- 
length]. These data were incomplete in the sense that 
they did not include the correct TT radiative piece in 
the metric, and thus did not contain realistic gravita- 
tional waves. 

In this paper, we revisit the PN data problem in 
ADM-TT coordinates, with the aim of supplying Numer- 
ical Relativity with initial BHB data that extend as far 
as necessary, and contain realistic gravitational waves. 
To do this, we have evaluated the "missing pieces" of 
Schafer's TT metric for the case of two non-spinning par- 
ticles. We have analyzed the near- and far-zone behavior 
of these data, and incorporated them numerically in the 
Cactus [4^ framework. In principle, the most accurate 
PN metric available could be used at this step, but it is 
not currently available in ADM-TT form. 

The remainder of this paper is laid out as follows. In 
Section |TT1 we summarize the results of Schafer (1985) 
[ill, and Jaranowski & Schafer (1997) [ill and their ap- 
plication by Tichy et al. (2003) [iH|, to the production 
of puncture data for numerical evolution. In Section 
mil we describe briefly the additional terms necessary 
to complete li^"^ to order (w/c)^, deferring details to the 
Appendix. In Section IIVI we study the full data both 
analytically and numerically. Section IVl summarizes our 
results, and lays the groundwork for numerical evolution 
of these data, to be presented in a subsequent article. 



precisely by a TT radiative part: 

9^j = (^1 + ^^) 
tt' = 0. 
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The fields (j), tt*-' and hj^^ can all be expanded in a post- 
Newtonian series. Solving the constraint equations of 
S-f-l general relativity in this gauge, [il, [i^ obtained ex- 
plicit expressions valid up to 0{v/cf' in the near zone, 
incorporating an arbitrary number of spinless point par- 
ticles, with arbitrary masses m^. For N particles, the 
lowest-order contribution to the conformal factor is^: 
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where = y/x — xa is the distance from the field point 
to the location of particle A. 



In principle hj^^ is computed from 
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where is the (flat space) inverse d'Alembertian (with 
a "no- incoming-radiation" condition [43|). s^i is a non- 
local source term and 5j^^^ is the TT-projection oper- 
ator. In order to compute hj^ we first rewrite Eq. (|4]) 
as 
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[A-i + (n-i-A-i)]sfc; 
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of h. 
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Note that the near-zone approximation h^^ „^ .^^^ 

has already been computed in [43t up to order 0(t;/c)^ 
(see also Eq. [TH below) . The last term in Eq. ^ is diffi- 
cult to compute because 
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is a non-local source. However, we can approximate Ski 
by 
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- 0{v/cf (8) 



II. ADM-TT GAUGE IN POST-NEWTONIAN 
DATA 



The "ADM-TT" gauge [43|, |47l is a 3-hl split of data 
where the three-metric differs from conformal flatness 



^ We explicitly include the gravitational constant G in all expres- 
sions here, as the standard convention G = 1 used in Numerical 
Relativity differs from the convention IGttG = 1 employed by 
[43,111. 
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in the near zone. Furthermore, outside the near zone The first term in (fTO|) . fij^^^'^^ can be expanded in v/c 
h^^/j- > ~ ^/r^, so that h'^^,,. . falls off much faster than as 

rest of ft-^"^, which falls off like 1/r. Hence 



h^-^ = h 



TT (NZ) rTT kl 
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where hj^^^^^^^^-^ can be neglected if we only keep terms up 

to 0(t;/c)^ generally, and 0{l/r) at infinity. 

The full expression for hj^ for N interacting point 
particles from Eq. (4.3) of |43| is: 
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The leading order term at 0(v/c)^, is given explicitly by 
(10) Eq. (A20) of Q: 



;,TT (4)., = ^ ^ _1_ I [II f _5 {fiA ■ PaY] S'' + 2p\p\ + [3{hA ■ pa)' - 5 II PA IP] n\ + 12(n^ • paW^Pa } 



4 ^ mATA 

A B^A 



-32 



rriAmB 



32 / 1 

SAB V^AS 



'AS 
7"^ / f A 



' AS 



nB \rB 



A"'AB ' 
17 



tabta 
4 



TABTA TATB 



1 

SAB 

1 1^^ 

""AB V^a 

-(- 

SAB \rA 

I — 



AS "AS 
2 



rA + J-s 



' AS 



12 

^AS 



3rA 

4 
rAB 







\rA 


SAB J _ 



'A "S 



'-A'M 



(12) 



where sab = ^a + 7's + tab- The other two terms 
in Eq. (jlOp can be shown to be small in the near zone 
(r ^ A, where the characteristic wavelength A ~ lOOA/ 
for TAs ~ 107\jf). However, ft,^"^'^^' is only a valid ap- 
proximation to hj!^ in the near zone, and becomes highly 
inaccurate when used further afield. 

Setting aside these far-field issues, Ticliy et al. [i^ ap- 
plied Schafer's formulation, in the context of a black-hole 
binary system, to construct initial data that are accurate 
up to 0{v/c)^ in the near zone. They noted that the 
ADM-TT decomposition was well-adapted to the use of 
a puncture approach to handle black-hole singularities. 
This approach is essentially an extension of the method 
introduced in [To| . It allows a simple numerical treatment 
of the black holes without the need for excision. 

The PN-based puncture data of Tichy et al. have 
not been used for numerical evolutions. This is in part 
because these data, just like standard puncture data [Tol . 
|49| . [sol , [sij , do not contain realistic gravitational waves 

in the far zone: h^^ does not even vaguely agree 
with the 2PN approximation to the waveform amplitude 



nor with the quadrupole approximation to the waveform 
phase for realistic inspiral. 

To illustrate this, we restrict to the case of two point 
sources, and compute the "plus" and "cross" polariza- 
tions of the near-zone approximation for hj^ : 

,(NZ) ,TT(NZ) i j 

,(NZ) ,TT(NZ) i j 

= ege'^. (14) 

For comparison, the corresponding polarizations of the 
quadrupole approximation for the gravitational-wave 
strain are given by (paraphrasing Eq. (3.4) of [52|): 

h+ = (l+cos^ 6l)(7rGAl/Gw)^''^cos($Gw), (15) 

r 

/i^ =l^cos0(7rGAl/Gw)'/'sin($Gw), (16) 
r 

where = v^^^ M is the "chirp mass" of the binary, 
given in terms of the total PN mass of the system M = 
nil + S'lid the symmetric mass ratio v = mim2/A/^. 
The angle is the "inclination angle of orbital angular 
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momentum to the line of sight toward the detector" ; that 
is, just the polar angle to the field point, when the binary 
moves in the x-y plane. <f>GW and /gw are the phase and 
frequency of the radiation at time t, exactly twice the 
orbital phase $(t — r) and orbital frequency ^{t — r) /2tt. 

The lowest-order PN prediction for radiation-reaction 
effects yields a simple inspiral of the binary over time, 
with orbital phasing given by 



V 

n(T) = ^_e-3/8, 

^ ' 8GM 



(17) 
(18) 



where = ly {tc ~ t)/5GM, M and ly are given below 
(fTHll , and tc is a nominal "coalescence time" . To evaluate 
(|13m4p . we need the transverse momentum p correspond- 
ing to the desired separation ri2. The simplest expression 
for this is the classical Keplerian relation, which we give 
parameterized by ^{t): 



6 = 7r/4, (f> = 0, for a binary in the x-y plane, with ini- 
tial separation ri2 = lOM. The orbital frequency of the 
binary is related to the separation ri2 and momenta p en- 
tering by p9ll20p . To this level of approximation, the 
binary has a nominal PN coalescence time tc ~ 78QM. 
As might have been anticipated, both phase and ampli- 

TT (4) 

tude of h^j are wrong outside the near zone. This 

TT (4) 

means that the data constructed from h-- have the 
wrong wave content, but nevertheless these data are still 
accurate up to order (v/c)^ in the far zone. 

It is evident from the present-time dependence of 
that it cannot actually contain any of the past history 
of an inspiralling binary. We would expect that a cor- 
rect "wave-like" contribution should depend rather on 
the retarded time of each contributing point source. It 
seems evident that the correct behavior must, in fact, be 
contained in the as-yet unevaluated parts of PU)) . The 
requisite evaluation is what we undertake in the next sec- 
tion. 



ri2 = Gi/3M(Mfl)-2/3, 
p = Mv(GAmfl^. 



(19) 
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In Fig. [T] we compare the plus polarization of the two 
waveforms (|13l) and ([T5|) at a field point r = 100 M, 



III. COMPLETING THE EVALUATION OF hj^ 

To move forward, we will simplify (jTU]) and to the 
case of only two particles. Then pU)) reduces to: 
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where 
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Here, the "TT projection" is effected using the operator 
:= 51 — ki y /k'^. For an arbitrary spatial vector u. 



U, + - 
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Details on the evaluarion^f ihese^terms are presentC' 
in Appendix [XI After calculation, we write the result as 



a sum of terms evaluated at the present field-point time 
t, the retarded time defined by 



i-t:,-r^(t:4) = o, 

and integrals between and t, 



(24) where the three parts are given by: 
ited 



(25) 



(26) 
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FIG. 1: Plus polarization of the quadrupole (black/solid) and 
near-zone (red/dashed) strains observed at field point r = 
fOOM, 6 = n /A, (fi — 0. The binary orbits in the x-y plane, 
with initial separation ri2 = lOM, and a nominal coalescence 
time tc ~ 780M. Both phase and amplitude of /i^^ are 
very wrong outside the near zone. 



H^ij, ^ [u] t\ t] 



1 G 



4 rA{t) 
+12 {u-nA)u^' 
G 



= -G dT 



t'; rA{Tf 



{ [-5 u^ + %{u- fiAf] S'^ +6u' - 12 {u ■ fiA) u'-' 



+ [9u^ - 15 {u-fiAf] n\n^^ 



G / dT 



I [li^ _ 5 . fiAf] 5« J +2u' - 20 (u • ua) u^' n^A 



[-5m^ + 35 (u • n^)^] n\n\^ 



(27) 



(28) 



(29) 



In Fig. [21 we show the retarded times calculated for 
each particle, as measured at points along the x axis, for 
the same orbit as in Fig. [T] Wc also show the corre- 
sponding retarded times for a binary in an exactly cir- 
cular orbit. Since the small-scale oscillatory effect of the 
finite orbital radius would be lost by the overall linear 
trend, we have multiplied by the orbital radius. 



A. Reconciling with Jaranowski & Schafer's h. 



From the derivation above it is clear that hJi includes 
retardation effects, so it will not depend solely on the 
present time. We might even expect that all "present- 
time" contributions should vanish individually, or should 
cancel out. It can be seen easily from (|27p that the "i" 
part of the second and third terms of Eq. (22) exactly 
cancel out the "kinetic" part (first line) of Eq. Thus, 
we can simply remove that line in Eq. (jl2[) . and use the 
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FIG. 2; Retarded times for particles 1 and 2, as measured by 
observers along the x axis at the initial time t = 0, for the 
binary of Fig. [1] To highlight the oscillatory effect of the 
finite-radius orbit on , we first divide by the average field 
distance r. 



also cancel the remaining, "potential" parts of Eq. (|12p . 
The answer is "not completely" ; expanding in powers of 
l/r, we find: 



-5{l + 6W^-7 W^) nu +2W{7 + 9W^) {n^, + nr2j nu)} + 0{l/r^)i31) 

I 



where W = sin 6* cos(0 — and is the orbital 

phase of particle 1 at the present time t. That is, the 
"new" contribution cancels the l/r and pieces of 

TT (A) 

h^j entirely. In the far zone the result is thus smaller 
than the lijj^^^^^^ term which we are ignoring everywhere, 
since it is small both in the near and the far zone [ist . 

We note here two general properties of the contribu- 
tions to the full hj"^ . 

'■J 

1. In the near zone hj^ ^^"^ is the dominant term since 

all other terms arise from — A^^)ski- Thus 

all other terms must cancel within the accuracy of 
the near-zone approximation. 

TT (A) 

2. h^j ^ ' is wrong far from the sources; thus, the 

TT (4:) 

new corrections should "cancel" h^j ^ ' entirely, far 
from sources. Note, however, that while hij = 
— □~g(S/ci depends only on retarded time, its TT- 
projection hj^^ = Sj^'^'^hki has a more complicated 
causal structure; E.g. the finite time integral comes 
from applying the TT-projection. [Proof: Even if 

we had a source given exactly by Ski , h^^ would 
depend only the present time, hij would depend 
only on retarded time, and hj/^ would (as we have 
computed) contain a finite time integral term.] 



Additionally, the full hj/^ agrees well with quadrupole 
predictions, which we demonstrate in Section [IVI 



IV. NUMERICAL RESULTS AND INVARIANTS 
A. Phasing and Post-Keplerian Relations 

It has been known for some time (see for example [ssj ) 
that gravitational wave phase plays an even more impor- 
tant part in source identification than does wave ampli- 
tude. In PN work, phase and amplitude are estimated 
somewhat separately; the amplitude requires knowledge 
of the time-dependent multipoles, used in developing the 
the full metric, while the phase can be relatively simply 
approximated from the orbital equations of motion, tak- 
ing into account the gravitational wave flux at infinity to 
evolve the orbital parameters jsj]. 

The quadrupole waveform introduced for the compar- 
ison in Fig. [T] had an amplitude accurate to 0{v/c)^ 
and the simplest available time evolution for the phase. 
Waveform phase is a direct consequence of orbital phase. 
To lowest order, we could have assumed a binary mov- 
ing in a circular orbit (of zero eccentricity) since, up to 
2PN order, we can have circular orbits, where the linear 
momentum, p, of each particle is related to the separa- 



tion ri2 by, say, Eq. (24) of [45[. Nevertheless, circular 
orbits are physically unrealistic - since radiation reac- 
tion will lead to inspiral and merger of the particles - 
and Eqs. (jlTHlSp already include leading-order radiation- 
reaction effects. Moreover, the phase errors that would 
accrue from using purely circular orbits would be larger, 
the further from the sources we tried to compute them. 

The calculations of section IIIII lead to waveform am- 
plitudes that are accurate at 0[v/c)^ everywhere. How- 
ever, we desire that our initial-data wave content already 
encode the phase as accurately as possible. Highly ac- 
curate phase for our initial data (via h'^'^), and hence in 
the leading edge of the waveforms we would extract from 
numerical evolution, is critical for parameter estimation 
following a detection. 

For demonstrative purposes, in this section, we will 
restrict ourselves to the simplest phasing relations con- 
sistent with radiation- reaction inspiral as given by Eqs. 
pTmsp . while using higher-order PN expressions than 
Eqs. pni-l^ for relating the orbit to the phase. For ex- 
ample, from [5^1, we have found to second PN (beyond 
leading) order: 



ri2(») 
GM 



Mv 
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(33) 



and we note that higher-order equivalents of these can 
be computed from [56| . 

In the numerical construction of initial data, the pri- 
mary input is the coordinate separation of the holes. In 
placing the punctures on the numerical grid, the separa- 
tion must be maintained exactly. To ensure this, we in- 
vert Eq. ((5^ to obtain the exact corresponding to our 
desired ri2. Then we use Eq. ^TE\\ with i = to find the 
coalescence time tc that yields this Qr- Once wc have ob- 
tained tc, we then find the orbital phase $ and frequency 

at any source time t directly from Eqs. (|17m8|) . and 
the corresponding separation ri2 and momentum p from 
Eqs. (|32ll33p . or their higher-order equivalents. 

In Fig. [21 we show a representative component of the 
retarded-time part of hj^"^ for both circular and leading- 
order inspiral orbits. For both orbits, we use the ex- 
tended Keplerian relations ([5^ and ([55)1 : otherwise the 
orbital configuration is that of Fig. [1] The coalescence 
time is now tc ^ llOOAf. Wc can see that the cumu- 
lative wavelength error of the circular-orbit assumption 
becomes very large at large distances from the sources. 
This demonstrates that using inspiral orbits instead of 
circular orbits will significantly enhance the phase accu- 
racy of the initial data, even though circular orbits are 
in principle sufficient when we include terms only up to 
0(u/c)^ as done in this work. From now on wc use only 
inspiral orbits. 




FIG. 3: The xx component of the full h^j for a binary with 
initial separation ri2 = 10 Af in a circular (black/solid) or 
inspiralling (red/dashed) orbit. Both fields have been rescaled 
by the observer radius r = z to compensate for the leading 
1/r fall-off. The orbital configuration is the same as for Fig. 
[Jl apart from the Keplerian relations, where we have used the 
higher-order relations (|32jl33|) . yielding tc ~ llOOAf. Note the 
frequency broadening at more distant field points. 




FIG. 4: Plus and cross polarizations of the strain observed at 
field point r = lOOAf, 6 ^ tt/4, = 0. Both the quadrupole- 
approximation waveform (black/solid and green/dot-dashed) 
and the full (red/dashed and blue/dotted) waveforms coming 
from hfj^ axe shown. The orbital configuration is the same as 
for Fig. \T\ 



Next, we compare our full waveform hjj^ (expressed 
as the combinations h+ and hx) at an intermediate-field 
position (r ~ lOOAf , 9 ~ 7r/4, = 0) to the lowest-order 
quadrupole result. In Fig. 31 the orbital configuration 
is the same as for Fig. [T] As one can see, both the 
-|- and X polarizations of our hj^^ agree very well with 
quadrupole results, as they should. We demonstrate the 
near- and intermediate-zone behavior of the new data on 
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FIG. 5: Plus and cross polarizations of the strain observed at 
t = along the z axis. We show the near-zone (solid/black), 
the quadrupole (dashed/red) and full (dot-dashed/green) 
waveforms. All waveforms have been rescaled by the observer 
radius r — z to compensate for the leading 1/r fall-off. The 
orbital configuration is the same as for Fig. [T] 



FIG. 6: Upper panel: Hamiltonian constraint violation along 
the y axis of our new data in the near zone, as a function of 
binary separation ri2. Lower panel: Momentum constraint 
(j/-component) violation of the same data along the x axis. 
The orbital configuration is that of Fig. [S] Distances have 
been scaled relative to ri2, so that the punctures are initially 
at y/ri2 = ±0.5. 



the initial spatial slice in Fig. [5] The quadrupole and 
full solutions agree very well outside ~ lOOAf . However, 
the full solution's phase and amplitude approach the NZ 
solution closer to the sources. 



B. Numerical Implementation 

After having confirmed that wc have a PN three-metric 
gij that is accurate up to errors of order 0{v/c)^, and 
that correctly approaches the quadrupole limit outside 
the near zone, we are now ready to construct initial data 
for numerical evolutions. In order to do so, we need the 
intrinsic curvature Kij, which can be computed as in 
Tichy et al. [i^ from the conjugate momentum. The 
difference is that here we use the full hj^ instead of the 



near-zone approximation h 
momentum 14 



TT (4) 



y to obtain the conjugate 
The result is 



-27 

^(3) 



2 



('^(2)7r 



~ij nTT 
(3)'' 



-0{v/c) 



(34) 



where the error term comes from neglecting terms like 



/jTT 

ij,{div) 



at 0{v/c) 



in hjj^ , and where ippN, t^^^^''^^ and 



An additional 



(/)(2) can be found in Tichy et al. [4 
difference is that the time derivative of hj,^ is evaluated 
numerically in this work. Note that the results for gij 
are accurate up to 0{v/c)'^, while the results for Kij are 



accurate up 0(v/c)^, because Kij contains an additional 
time derivative [H, [13, [Ell . 

Next we show the violations of the Hamiltonian and 
momentum constraints computed from gij and Kij, as 
functions of the binary separation ri2. As we can see in 
both panels of Fig. El the constraints become smaller for 
larger separations, because the post-Newtonian approxi- 
mation gets better. Note that, as in [i^, the constraint 
violation remains finite everywhere, and is largest near 
each black hole. 



C. Curvature Invariants and Asymptotic Flatness 

In analysis of both initial and evolved data, it is often 
instructive to investigate the behavior of scalar curva- 
ture invariants, as these give some idea of the far-field 
properties of our solution. We expect, for an asymptoti- 
cally fiat space-time, that in the far field, the speciality 
index S = TlJ^jX^ will be close to unity. This can 
be seen from the following arguments. Let us choose a 
tetrad such that the Weyl tensor components i/ji and ^/'a 
are both zero. Further, we assume that in the far field 
'00 £md ■04 are both perturbations of order e off a Kerr 
background. Then 



1-3 



0004 
01 



(35) 
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which is indeed close to one. Note however, that this 
argument only works if the components of the Weyl ten- 
sor obey the peeling theorem, such that ip2 ^ 0{r~^), 
ipo ^ 0{r~^) and ip4 ~ 0{r~^). In particular, if ipo falls 
off more slowly than 0(r~^), S will grow for large r. 
Now observe that ipo ^ 0{r^^) ^ M^/r^ is formally of 
0{v/c)^. Thus, in order to see the expected behavior of 
iS ~ 1 in the far- field we need to go to 0{v/c)^. If wc only 
go to 0(u/c)^ (as done in this work) ipo consists of un- 
controlled remainders only, which should in principle be 
dropped. When we numerically compute S we find that 
for our data, S deviates further and further from unity 
for large distances from the binary. This reflects the fact 
that the so-called "incoming" Weyl scalar 2po only falls 
off as due to imcontroUed remainders at 0{v/c)^, 

which arise from a mixing of the background with the 
TT waveform. 



nation of the gravitational radiative modes; it may be 
possible to produce them in ADM-TT coordinates via 
a contact transformation, or by direct calculation (see, 
e.g. llOl)- For initial separations similar to the fiducial 
test case of this paper, ri2 = lOM, the order necessary 
for clean matching of the initial wave content with the 
new radiation generated in evolution should not be par- 
ticularly high j26|. As noted, the Kcplcrian relations 
Eqs. (|32II33[) can easily be extended to higher PN order. 

The data presented already allow for arbitrary initial 
mass ratios this introduces the possibility of significant 
gravitational radiation in odd-/ multipoles, together with 
associated phenomena, such as in-plane recoil "kicks" . 
An interesting future development of these data will be 
the inclusion of spin angular momenta on the pre-merger 
holes. This will open our initial-data prescription to de- 
scribing an even richer spectrum of binary radiation. 



V. DISCUSSION AND FUTURE WORK 



Exploring and validating PN inspiral waveforms is cru- 
cially important for gravitational-wave detection and for 
our theoretical understanding of black- hole binaries. Our 
goal has been to provide a step forward in this under- 
standing by building a direct interface between the PN 
approach and numerical evolution, along the lines ini- 
tially outlined in Ref. [i^. In this paper we have 
essentially completed the calculation of the transverse- 
traceless part of the ADM-TT metric to 0(w/c)* pro- 
vided in [45, yielding data that, on the initial Cauchy 
slice, will describe the space-time into the far-field. 
We have incorporated this formulation into a numerical 
initial-data routine adapted to the "puncture" topology 
that has been so successful recently, and have explored 
these data's numerical properties on the initial slice. 

Our next step is to evolve these data with moving 
punctures, and investigate how the explicit incorporation 
of post-Newtonian waveforms in the initial data affects 
both the ensuing slow binary inspiral of the sources and 
the release of radiation from the system. We note es- 
pecially that our data are non-confornially flat beyond 
0{v/c)^. We expect our data to incorporate smaller 
unphysical initial distortions in the black holes than is 
possible with conformal flatness, and hence less spurious 
gravitational radiation during the numerical evolution. 
We see this as a very positive step toward providing fur- 
ther validation of numerical relativity results for multiple 
orbit simulations, since it permits comparison with PN 
results where they are expected to be reliable. Our initial 
data will also allow us to fully evaluate the validity of PN 
results for merging binaries by enabling comparison with 
the most accurate numerical relativity results. 

We expect that further development of these data will 
certainly involve the use of more accurate orbital phas- 
ing information than the leading order given by Eqs. (|17t - 
llSp . This information is available in radiative coordinates 
(see, e.g. Eq. (6.29) of [5^) appropriate for far-field eval- 
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APPENDIX A: DETAILS OF INTEGRAL 
CALCULATION 



Here we present some more details of the calculations 
that lead to the three contributions to Eq. ([^5]) : Eqs. 
(|27ll29p . Inserting Eq. ^ in the general integral 



we can write H^j,^[u\ as a combination of scalar and 
tensor terms: 



H. 



TTAr 



16ttG I dr 
,2 



Ia 



+ 



lij A + 



UcUd 



red r 
^A "ij 



UcUd 
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where the "/" integrals are defined as: 



Ia 



— 



^A 



(Pkdw {uj/k) 



2 Ak TA COS 9—i uj T 



(2 7r)4 k^-{uj + ief 

d^kduo k k' 
(2 7r)4 ~fc2~ 

(cl'/A;)'^ gifc^A cos 9— iu)T 
fc2 - (cj + ie)2 ' 
d^kduok k^ k'= k'^ 

k'^ — {uj + i e)2 



(A2) 



(A3) 



(A4) 



Here T — t — t, and Fyi = a? — xa- We have also taken 
our integration coordinates such that ta lies in the z di- 
rection, so that the dummy momentum vector k satisfies 



So the next integrals will differ in their 9 dependence, 
contained in the powers of w above. The 9 integrals will 
contain the following basic types: 



92{a) 
54(a) 



dwe" 



sinha 



(A7) 



-1 
+1 



+4- 



sinh ( 



sinh a cosh a 
2 4 — 



(A8) 



+1 



dw w e 



4 a w 



sinh a cosh a 
2 8 TT- 



sinh a cosh a 
-24 48 :r- 



sinha 



(A9) 



Now and can be written as the linear combi- 

nations: 



k ■ -Fa = kvA cos 9, (A5) 
d^k = k^dk s\n9d9d(t). (A6) 

Define the unit orthogonal vectors (0,0,1) , 1 = 

(cos^, sin0, 0). Then we can write 

k = k cos 9 fiA + k sin ^ =^ k ■ va ~ ^Ak ■ fiA ■ 
We can also define a projector tensor onto I: 

1. Angular integration 

We will neglect the A subscript for now, until it be- 
comes relevant again. To calculate the integrals (jA2IIA4p , 
we begin with the (f> integration. The only (j) dependence 
comes from the £ parts of the k terms. It can be seen 
from elementary trigonometric integrals that: 



^r^^r = 0, j dipt = t:Q''^, 
jdiljei^tl"^ = J (Q°''g"'^-f Q'''=Q'"* + Q"''Q''") . 



We use these to calculate the integrals for Ij^ and 



I'X^'^'K Define w = cos 9. Then 



'1 = 2 7r, 



k'' k^ 
k" k^ ¥ k'^ 

¥ 



2TTw'^n''n'' + tt{1- w^)Q 



2\ r^ab 



+6 7ru;2(i„u;2)Q(afc„c^d 

+ ^(l_y,2)2g(a6Qcd)^ 



ja b 
ja bed 



1 



(AlO) 



8 

L J T 

/ here can be expressed in terms of go (a) above: 

du r d^k {uj/kf .krcosB-ru:T 



(All) 



(2 7r)3 J 2 7r fc2 - (cj + ie)2 



d(jj 
(2^ 

duj 
(2^ 



,2 „-iujT 



dk 



1/2 



fc2 ~ (uj + i e)' 



Jo 



■goiikr) 
(A12) 



The 1/2 factor is because we moved to integrating k over 
the whole real line instead of the positive half-line (this 
is permissible as gn(a) is an even function of a). K and L 
are defined analogously to /, but with extra even powers 
of cos 9 ~ w: 



K 



duj f d^k {uj/kf 



(2 7r)3 J 2 7r fc2 - (cj + ie)- 



k r cos 8 — iu T j 



duo 

(2^ 

duj 



,2 „-iwT 



dk 



1/2 



k^ — (lu + i e) 



2 92(1 kr) 



J2 



L = 



(2 7r)3 

duj f d^k {oj/kf 



(A13) 



(2 7r)3 J 2 7r fc2 - (cJ + ^e)■ 



, gifercose-i<^T^,Qg4| 



duj 
(2^ 

dio 
(2^ 



2 -iiuT 

UJ e 



, ,2 -ILOT J 



dk 



1/2 



k'^ ^ (uj + i e)2 



34 (« kr) 



(A14) 
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2. Momentum integration 



Now we address the k integrals, defined as: 



dkf+{k)+ dkf-{k), 



dk fn{k) 



where we collect the positive exponents in the gn in the 
integrand of f^{k), and the negative exponents in f^{k): 



, fn{k) 



9niikr)/2 
k^ - {uj + iey 



We calculate this as the sum of contour integrals of the 
"plus" and "minus" integrands (necessary, as the oppo- 
site signs require different contours). Each of these has 
poles at /c = 0, k ~ k-^- = uj + i e, and k = k^ = —lu — ie 
(the first of these is from the gn)- We integrate the "plus" 
integrands anticlockwise around the contour Ci , and the 
"minus" integrands anticlockwise around the contour C2 
(see Fig. [7]); taking the limit |fc| 00, the contribu- 
tion from the curved segments vanishes, and the residue 
theorem gives us: 

Jn = 2TTiRes[fn ,k+] — 2niRes[fn ,k^] 

+TriRes[f+,0]-niRes[f-,0]. (A15) 



Calculating the residues, we find the values of each of 
the Jn- 



Jo 



r (uj+i e) 



(A16) 



^^tr(u>+te) ^^ir{^+ie) [_2 + 2 i r (w -|- i e)] 



r {uj + i tY 
2tx 

-3 r'^ {uj + i (if + ir^ {uj + i ef 

24 TT 



(A17) 

[6 — 6 i r (lj -f i e) 

(A18) 



3. Frequency integration 

Now we perform the uj integration. Inserting the re- 
sults (|A16IIA18|) into (|A12IIA14p respectively, we see that 
each of /, K and L contains a delta function, which we 



can extract: 
1 

/ 
K 



4 TT r 
1 

4 TT r 



[5{T-r)-5{T)l 



6(T -r) + e 



dio 



L = 



4 TT r 



dhJ 

J2nf^ 

5{T - r) 
duj 



(2 7r)3 
doj 



{2nf 



(2 7r)3 

where the new terms on the right-hand side come from 
the Jn above, grouped by exponential, as that is what 
determines the contours chosen during integration (see 
Fig. ED: 



F2a{uj) 
F2b{uj) 
Fiai^) 



TT uj"^ [~2 + 2ir [u + i e) 

2'KUj'^ 



r3 {uj + 16)-^ 



[24- 24zr(u) + ie) 



(tj + i e)^ 

-12 (w + je)^ + 4ir^ (w + ie)^] , 
247rtj2 



F^b{uj) = - , ■ \e,- 

r° {uj ^ I e)° 

Now the residues are as follows (taking the e ^ limit): 

27riT 



Res 



Res 



e-'"(^-'^)F2,(c.),-ze 
Res [e-'"^F2fc(w),-ze 
-'"(^-'■)F4a(c.),-ze' 



Res [e-'"^ 



F4f,(w), -ie] 



2'KiT 
47riT3 



The only pole is at w = —i e, so if we can close the contour 
in the upper half-plane, we will get zero. 

• For T < 0, both the "a" and "b" integrals can be 
closed in Ci. Result: zero contribution. 

• For < T < r, the "a" integrals can be closed 
in Ci, but the "b" integrals must be closed in C2. 
Result: "b" contribution. 

• For T > r, both the "a" and "b" integrals must be 
closed in C2. But then the "a" and "b" residues 
cancel out. Result: zero contribution. 

Thus the only interesting contribution happens in the 
interval 0<T<r<^t — r(T) < t < t. In this case, the 
final integrals yield 



duj 
(2^ 

duj 
(2^ 



-iuj iT-r) 



F2b{i0) 



T 



2 7rr3 ' 



rp3 

-iuj (T-r) p ,^ _ 

' rib[i^) — 



12 





^\ y 

u + ze\ 
• 










( 

1 V. 
\ • 

\ —u — ie 


> 


k 






UJ 



FIG. 7: Contours needed to complete integration over k (left) and uj (right). 

leading to the final result for K and L: We use these to calculate the and 

I = J_J(T-r)--^<5(r), 
47rr 47rr 

K = ^5{T-t)-Q{T)Q{t-T)-^, 
4 TT r TTr° 



jijkl 



5{T- 



<d{T)Q{r-T) 



T 



2 7: 



1 



^ 5(T) + e(r)e(r-r) ^ 



n'^ n' 



^—5{T-r)-Q{T)Q{r-T)^ 



2 \ 47rr 



2 TT r"^ 
T T3 



(M9) 



1 f T 

-^5{T) + Q{T) e(r - T) — - — 



(A20) 



4. Time integration 



The final integrations will be over the source time t. The "crossing times" for the two functions are t ^ t and 
T ~ f , where t is the present field time, and V the corresponding retarded time defined by (j25p . Now taking a general 
function ^(t). wc find that 



/oo 
dTlAvir) = 
-00 

drrjy{r) = 

00 



y{t) 



v{t\) 



"f^A "IT- A Z 

^ 4 TT r A 



Ua n\ n\ riA ^ 

4 TT TA 
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4:77 rA^r)^ 



3 n(^j n^i) vi.T) 



^Qa Qa 4^^^ 
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These can now be substituted into the general integral (|Aip . We write the result as a sum of terms at the present 
field-point time t, the retarded time i^, and interval terms between them, 

Hjip ^ [li] = i7rp-?p ^ [u; t\ + _ffrp^ ^ [u; i^] + H^r^ ^ [u; ^ , 



^TT A I"*^' 
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+ 



t 




d'^ 

2 



n'X n\ S"' 



n\ n\ n\ n\ 



Hrflj, ^ [U] t^ 



t] = -AG f dr^lj^ 



o i j 
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